Short Time Fourier Transform (STFT), Wigner Distribution, and Pseudo Wigner Distribution

last update: 9/30 (2017)



In [1]:

    
VERSION









    Out[1]:





v"0.6.1"



In [2]:

    
import DSP
using PyPlot



In [3]:

    
# include all modules in juwvid
include("../juwvid.jl")









    Out[3]:





juwvid

Single Component Linear Signal



In [4]:

    
# Figure 1
nsample=1024
x,y,iw,ynorm=sampledata.genlinfm(nsample,1.0,0.01);
fig=PyPlot.figure(figsize=(10,3))
PyPlot.plot(x,y)
PyPlot.savefig("linearFMnoise.png")

Short Time Fourier Transform



In [5]:

    
tfrst=stft.tfrstft(y);









    



Use fft.



In [6]:

    
# Figure 2
a=juwplot.wtfrshow(abs.(tfrst).^2,x[2]-x[1],x[1],x[end],NaN,NaN,0.7,"CMRmap")
PyPlot.xlabel("time")
PyPlot.ylabel("frequency")
PyPlot.title("Spectrogram")









    












    Out[6]:





PyObject Text(0.5,1,'Spectrogram')

Changing the window size



In [7]:

    
tfrst4=stft.tfrstft(y,NaN,NaN,NaN,NaN,NaN,4);
tfrst8=stft.tfrstft(y,NaN,NaN,NaN,NaN,NaN,8);
tfrst16=stft.tfrstft(y,NaN,NaN,NaN,NaN,NaN,16);
tfrst32=stft.tfrstft(y,NaN,NaN,NaN,NaN,NaN,32);









    



Use fft.
Use fft.
Use fft.
Use fft.



In [8]:

    
# Figure 3
fig=PyPlot.figure()
ax = fig[:add_subplot](2,2,1)
a=juwplot.wtfrshow(abs.(tfrst4).^2,x[2]-x[1],x[1],x[end],NaN,NaN,0.7,"CMRmap")
PyPlot.ylabel("frequency")
PyPlot.title("width = N/4")
ax = fig[:add_subplot](2,2,2)
a=juwplot.wtfrshow(abs.(tfrst8).^2,x[2]-x[1],x[1],x[end],NaN,NaN,0.7,"CMRmap")
PyPlot.title("width = N/8")
ax = fig[:add_subplot](2,2,3)
a=juwplot.wtfrshow(abs.(tfrst16).^2,x[2]-x[1],x[1],x[end],NaN,NaN,0.7,"CMRmap")
PyPlot.xlabel("time")
PyPlot.ylabel("frequency")
PyPlot.title("width = N/16")
ax = fig[:add_subplot](2,2,4)
a=juwplot.wtfrshow(abs.(tfrst32).^2,x[2]-x[1],x[1],x[end],NaN,NaN,0.7,"CMRmap")
PyPlot.xlabel("time")
PyPlot.title("width = N/32")









    












    Out[8]:





PyObject Text(0.5,1,'width = N/32')

Comparison with the (pseudo) Winger Ville distribution



In [9]:

    
### generating the analytic signal of y and computing the WV and the pseudo WV
zs=DSP.Util.hilbert(y); 
tfrs=cohenclass.tfrwv(zs);
tfrps=cohenclass.tfrpwv(zs);









    



Single Wigner Ville
Use fft.
Single pseudo Wigner Ville
Use fft.



In [10]:

    
# Figure 4
fig=PyPlot.figure()
ax = fig[:add_subplot](2,2,1)
a=juwplot.wtfrshow(abs.(tfrst).^2,x[2]-x[1],x[1],x[end],NaN,NaN,0.7,"CMRmap")
PyPlot.ylabel("frequency")
PyPlot.title("Spectrogram")
ax = fig[:add_subplot](2,2,2)
a=juwplot.tfrshow(real(tfrs),x[2]-x[1],x[1],x[end],NaN,NaN,0.7,"CMRmap")
PyPlot.title("Wigner Ville")
ax = fig[:add_subplot](2,2,3)
a=juwplot.tfrshow(real(tfrps),x[2]-x[1],x[1],x[end],NaN,NaN,0.7,"CMRmap")
PyPlot.ylabel("frequency")
PyPlot.xlabel("time")
PyPlot.title("Pseudo Wigner Ville")
ax = fig[:add_subplot](2,2,4)
PyPlot.plot(x,y,color="gray",lw=0.3)
PyPlot.xlim(0,x[end])
PyPlot.ylim(-1.2,1.2)
PyPlot.xlabel("time")
PyPlot.title("Input data")









    












    Out[10]:





PyObject Text(0.5,1,'Input data')

with noise



In [11]:

    
# Figure 5
using Distributions
d = Normal()
ynoise=y+rand(d,nsample)*std(y)*1.5
fig=PyPlot.figure(figsize=(10,3))
PyPlot.plot(x,ynoise)
PyPlot.savefig("linearFMnoise.png")



In [12]:

    
tfrst=stft.tfrstft(ynoise);
zs=DSP.Util.hilbert(ynoise); 
tfrs=cohenclass.tfrwv(zs);
tfrps=cohenclass.tfrpwv(zs);









    



Use fft.
Single Wigner Ville
Use fft.
Single pseudo Wigner Ville
Use fft.



In [13]:

    
# Figure 6
fig=PyPlot.figure()
ax = fig[:add_subplot](2,2,1)
a=juwplot.wtfrshow(abs.(tfrst).^2,x[2]-x[1],x[1],x[end],NaN,NaN,0.7,"CMRmap")
PyPlot.ylabel("frequency")
PyPlot.title("Spectrogram")
ax = fig[:add_subplot](2,2,2)
a=juwplot.tfrshow(abs.(tfrs),x[2]-x[1],x[1],x[end],NaN,NaN,0.7,"CMRmap")
PyPlot.title("Wigner Ville")
ax = fig[:add_subplot](2,2,3)
a=juwplot.tfrshow(abs.(tfrps),x[2]-x[1],x[1],x[end],NaN,NaN,0.7,"CMRmap")
PyPlot.ylabel("frequency")
PyPlot.xlabel("time")
PyPlot.title("Pseudo Wigner Ville")
ax = fig[:add_subplot](2,2,4)
PyPlot.plot(x,ynoise,color="gray",lw=0.3)
PyPlot.xlim(0,x[end])
PyPlot.ylim(-3.5,3.5)
PyPlot.xlabel("time")
PyPlot.title("Input data")









    












    Out[13]:





PyObject Text(0.5,1,'Input data')

Single component nonlinear signal



In [14]:

    
# Figure 7
xs,ys,iws,ynorms=sampledata.genfm(nsample,2*pi,2*pi/365.0,100.0,365.0);
fig=PyPlot.figure(figsize=(10,3))
PyPlot.plot(xs,ys)
PyPlot.savefig("nonlin.png")



In [15]:

    
tfrst=stft.tfrstft(ys);
zs=DSP.Util.hilbert(ys); 
tfrs=cohenclass.tfrwv(zs);
tfrps=cohenclass.tfrpwv(zs);









    



Use fft.
Single Wigner Ville
Use fft.
Single pseudo Wigner Ville
Use fft.



In [16]:

    
# Figure 8
fig=PyPlot.figure()
ax = fig[:add_subplot](2,2,1)
a=juwplot.wtfrshow(abs.(tfrst).^2,x[2]-x[1],x[1],x[end],NaN,NaN,0.7,"CMRmap")
PyPlot.ylabel("frequency")
PyPlot.title("Spectrogram")
ax = fig[:add_subplot](2,2,2)
a=juwplot.tfrshow(abs.(tfrs),x[2]-x[1],x[1],x[end],NaN,NaN,0.7,"CMRmap")
PyPlot.title("Wigner Ville")
ax = fig[:add_subplot](2,2,3)
a=juwplot.tfrshow(abs.(tfrps),x[2]-x[1],x[1],x[end],NaN,NaN,0.7,"CMRmap")
PyPlot.ylabel("frequency")
PyPlot.xlabel("time")
PyPlot.title("Pseudo Wigner Ville")
ax = fig[:add_subplot](2,2,4)
PyPlot.plot(xs,ys,color="gray",lw=0.3)
PyPlot.xlim(0,x[end])
PyPlot.ylim(-1.2,1.2)
PyPlot.xlabel("time")
PyPlot.title("Input data")









    












    Out[16]:





PyObject Text(0.5,1,'Input data')

multicomponent nonlinear signal



In [17]:

    
# Figure 9
nsample=512
t,ym=sampledata.genmultifm623(nsample);
fig=PyPlot.figure(figsize=(10,3))
PyPlot.plot(t,ym)
PyPlot.savefig("multiFM.png")



In [18]:

    
tfrst=stft.tfrstft(ym);









    



Use fft.



In [19]:

    
### generating the analytic signal of y
zm=DSP.Util.hilbert(ym); 
tfrs=cohenclass.tfrwv(zm);
tfrps=cohenclass.tfrpwv(zm);









    



Single Wigner Ville
Use fft.
Single pseudo Wigner Ville
Use fft.



In [20]:

    
# Figure 10
fig=PyPlot.figure()
ax = fig[:add_subplot](2,2,1)
a=juwplot.wtfrshow(abs.(tfrst).^2,t[2]-t[1],t[1],t[end],NaN,NaN,0.7,"CMRmap")
PyPlot.ylabel("frequency")
PyPlot.title("Spectrogram")
ax = fig[:add_subplot](2,2,2)
a=juwplot.tfrshow(abs.(tfrs),t[2]-t[1],t[1],t[end],NaN,NaN,0.7,"CMRmap")
PyPlot.title("Wigner Ville")
ax = fig[:add_subplot](2,2,3)
a=juwplot.tfrshow(abs.(tfrps),t[2]-t[1],t[1],t[end],NaN,NaN,0.7,"CMRmap")
PyPlot.xlabel("time")
PyPlot.title("Pseudo Wigner Ville")
ax = fig[:add_subplot](2,2,4)
PyPlot.plot(t,ym,color="gray",lw=0.3)
PyPlot.xlim(t[1],t[end])
PyPlot.ylim(-2.2,2.2)
PyPlot.xlabel("time")
PyPlot.title("Input data")









    












    Out[20]:





PyObject Text(0.5,1,'Input data')

For a such multimode signal, you can use S-method or higher order Wigner distribution instead of PWV.



In [ ]: